classdef State
    % State Class

    % The "state" or "fundamental" is the (random) probability of success 
    % of the bank's risky project. In the numerical examples, we assume that
    % the logit of the state has a normal distribution.

    properties
        Mean   = 0.25      % Mean of normal distribution
        Vol    = 0.50      % Standard deviation of normal distribution
    end

    methods (Static)
        function state = getState(logit)
            % The logistic function maps the logit into the state
            state = 1./(1+exp(-logit));
        end

        function [logit,Dlogit] = getLogit(state)
            % The "logit" maps the state into a normally distributed r. v.
            logit    = log(state./(1-state));
            Dlogit   = 1./(state.*(1-state));
        end
    end
    
    methods
        function CDF = cdfLogit(obj,logit)
            % CDF of the "logit": Y = log(state/(1-state))
            CDF = normcdf(logit,obj.Mean,obj.Vol);
        end

        function pdf = pdfLogit(obj,logit)
            % pdf of the "logit": Y = log(state/(1-state))
            pdf = normpdf(logit,obj.Mean,obj.Vol);
        end

        function CDF = cdfProb(obj,state)
            % CDF of the state
            logit   = obj.getLogit(state);
            CDF     = obj.cdfLogit(logit);
        end

        function pdf = pdfProb(obj,state)
            % pdf of the state
            [logit,Dlogit] = obj.getLogit(state);
            pdf = obj.pdfLogit(logit).*Dlogit;
        end

        function [nod,wgt] = quadNormal(obj,varargin)
            % This function performs a numerical quadrature of a normal r. v. 

            p = inputParser;
            p.addParameter('BoundLeft',[],@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('BoundRight',[],@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('CutOffs',[-7,7],@(x)isnumeric(x)&&length(x)==2);
            p.addParameter('NumPoints',2000,@(x)isnumeric(x)&&isscalar(x));
            p.parse(varargin{:});
            
            mu      = obj.Mean;
            sig     = obj.Vol;
            a       = p.Results.BoundLeft;
            b       = p.Results.BoundRight;
            m       = p.Results.CutOffs;
            N       = p.Results.NumPoints;
            
            amin = mu+m(1)*sig; bmax = mu+m(2)*sig;
            if isempty(a); a=amin; end
            if isempty(b); b=bmax; end

            cdf = @(x)cdfLogit(obj,x);
            grd = linspace(max(a,amin),min(b,bmax),N+1)';
            nod = (grd(1:N)+grd(2:N+1))/2;
            tmp = cdf(grd);
            wgt = (tmp(2:N+1)-tmp(1:N));
            if (a>amin)
                if (b<bmax)
                    wgt = wgt/sum(wgt)*(cdf(b)-cdf(a));
                else
                    wgt = wgt/sum(wgt)*(1-cdf(a));
                end
            else
                if (b<bmax)
                    wgt = wgt/sum(wgt)*cdf(b);
                else
                    wgt = wgt/sum(wgt);
                end
            end
        end

        function [nod,wgt] = quadState(obj,varargin)
            % This function performs a numerical quadrature of the state
            
            p = inputParser;
            p.addParameter('BoundLeft',[],@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('BoundRight',[],@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('CutOffs',[-7,7],@(x)isnumeric(x)&&length(x)==2);
            p.addParameter('NumPoints',2000,@(x)isnumeric(x)&&isscalar(x));
            p.parse(varargin{:});
            
            p_a = p.Results.BoundLeft;
            p_b = p.Results.BoundRight;
            m   = p.Results.CutOffs;
            N   = p.Results.NumPoints;
            mu  = obj.Mean;
            sig = obj.Vol;
            
            if isempty(p_a); a=mu+m(1)*sig; else; a=obj.getLogit(p_a); end                
            if isempty(p_b); b=mu+m(2)*sig; else; b=obj.getLogit(p_b); end                
            
            [logit_nod,wgt] = obj.quadNormal(...
                'BoundLeft',a,'BoundRight',b,'CutOffs',m,'NumPoints',N);
            nod = obj.getState(logit_nod);
        end        
    end    
    
end